# SI Figure 3
rm(list=ls())

# Load RCPs
## First column is temperature time series from deterministic EBM. initialised at 1C in 2020. Next 10,000 columns are temperature time series from stochastic EBM. (See Matlab script.)
# path <- "~/Dropbox/StochasticIAMS/"
path <- "~/Documents/Projects/IAM_Stochastic/"
files <- list.files(path=path, pattern="output.txt")
files <- c(files,list.files(path=path, pattern="Sigma.txt"))
d <- data.frame()
for (i in 1:length(files)){
	ds <- read.table(paste(path, files[i], sep=""), sep=",", header=FALSE)
	ds$rcp <- files[i]
	d <- rbind(d, ds)
}

# RGB colors for RCPs from IPCC AR5.
d$r <- NA
d$r[grep("26",d$rcp)] <- 0
d$r[grep("45",d$rcp)] <- 121
d$r[grep("60",d$rcp)] <- 255
d$r[grep("85",d$rcp)] <- 255
d$g <- NA
d$g[grep("26",d$rcp)] <- 0
d$g[grep("45",d$rcp)] <- 188
d$g[grep("60",d$rcp)] <- 130
d$g[grep("85",d$rcp)] <- 0
d$b <- NA
d$b[grep("26",d$rcp)] <- 255
d$b[grep("45",d$rcp)] <- 255
d$b[grep("60",d$rcp)] <- 45
d$b[grep("85",d$rcp)] <- 0

d$sigmaQ <- NA
d$sigmaQ[grep("lowSigma",d$rcp)] <- "low"
d$sigmaQ[grep("highSigma",d$rcp)] <- "high"
d$sigmaQ[grep("output.txt",d$rcp)] <- "mid"

# Ensemble size
nn = ncol(d) - 7
names(d) <- c("year","det",paste("X",seq(1:nn),sep=""), "rcp","r","g","b","sigmaQ")
d$year <- d$year + 2019

# Length of time series
n=length(unique(d$year))

# Economic assumptions 
gdp = 80*10^12 # initial aggregate consumption (80 trillion US dollars)
g = 0.019 # annual growth rate of undisturbed aggregate consumption
eta = c(0,1.45) # marginal utility of consumption
rho = c(0.04255,0.015) # pure rate of time preference (1.5% in DICE-2016)

s <- 1 # Discounting assumptions: r = 4.255 + 0 * 1.9 = 4.255

# Demographic assumptions
pop = 7.5 # billion people
pop.asympt = 11.5 # Asymptotic population, in billion
param = 0.134/5 # DICE-2016 "Population 2050 parameter"
pop.ts <- rep(NA,n)
pop.ts[1] <- pop
for (i in 2:n) { pop.ts[i] <- pop.ts[i-1]*(pop.asympt/pop.ts[i-1])^param}
pop.ts <- matrix(rep(pop.ts,nn+1),ncol=nn+1,nrow=n)
pop.ts <- pop.ts*10^9 # Population

# Undisturbed output per capita
gdp.ts <- gdp*(1+g)^seq(0,n-1,1)
gdp.ts <- matrix(rep(gdp.ts,nn+1),ncol=nn+1,nrow=n)
gdp.ts <- (gdp.ts/pop.ts)/1000 # undisturbed consumption per capita expressed in thousands of dollars.

# Discount factors
df.2 <- 1/((1+rho[s])^seq(0,n-1,1))
df.2 <- matrix(rep(df.2,nn+1),ncol=nn+1,nrow=n)

x <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp","r","g","b","sigmaQ"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.2 * pop.ts)
	x <- rbind(x, x.Weitzman)
}
x <- t(x)
colnames(x) <- files

## Scaling factors for plotting
UofC <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar <- UofC-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.2[,1]) # NPV of utility of undisturbed aggregate consumption

# Estimate risk premium
s <- 2 # Discounting assumptions: r = 1.5 + 1.45 * 1.9 = 4.255

# Discount factors
df.1 <- 1/((1+rho[s])^seq(0,n-1,1))
df.1 <- matrix(rep(df.1,nn+1),ncol=nn+1,nrow=n)

x1 <- data.frame()
for (j in files) {
	t <- d[d$rcp==j,]
	t <- t[,!(names(t) %in% c("year","rcp","r","g","b","sigmaQ"))]
	# Isoelastic utility
	x.Weitzman <- colSums((((gdp.ts*(1/(1+((t)/20.46)^2 + ((t)/6.081)^6)))^(1-eta[s]))/(1-eta[s])) * df.1 * pop.ts)
	x1 <- rbind(x1, x.Weitzman)
}
x1 <- t(x1)
colnames(x1) <- files

## Scaling factors
UofC1 <- ((gdp.ts[1,1])^(1-eta[s]))/(1-eta[s]) # Utility of undisturbed per capita consumption today (in thousands of dollars)
scalar1 <- UofC1-(((gdp.ts[1,1]-0.001)^(1-eta[s]))/(1-eta[s])) # Utility value of the marginal dollar
UofC1 <- sum((((gdp.ts[,1])^(1-eta[s]))/(1-eta[s]))*pop.ts[,1]*df.1[,1]) # NPV of utility of undisturbed aggregate consumption

# Plot
quartz(width = 10, height = 5, type = "pdf", "SI_figure_sigma", file = "SI_figure_sigma.pdf")
par (fig=c(0,1,0,1), # Figure region in the device display region (x1,x2,y1,y2)
	omi=c(.5,.7,0,0), # global margins in inches (bottom, left, top, right)
	mai=c(0,.1,0.4,.1)) # subplot margins in inches (bottom, left, top, right)
	layout(matrix(c(1:12), 3, 4, byrow = FALSE), widths=c(3,3,3,3), heights=c(3,3,3))

# Plot damage distributions
j = "RCP26outputhighSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="RCP2.6",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(10*10^12,60*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*60*10^12, y=0.9*max(dens$y), labels="Risk Premium", cex=1.2, pos=2)
text(x=0.99*60*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*60*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP26output.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(10*10^12,60*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*60*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*60*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP26outputlowSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(10*10^12,60*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*60*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*60*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

# Plot damage distributions
j = "RCP45outputhighSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="RCP4.5",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(30*10^12,120*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*120*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*120*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP45output.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(30*10^12,120*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*120*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*120*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP45outputlowSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(30*10^12,120*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*120*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*120*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

# Plot damage distributions
j = "RCP60outputhighSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="RCP6.0",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(50*10^12,220*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*220*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*220*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP60output.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(50*10^12,220*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*220*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*220*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP60outputlowSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(50*10^12,220*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*220*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*220*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

# Plot damage distributions
j = "RCP85outputhighSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="RCP8.5",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(300*10^12,680*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*680*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*680*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP85output.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(300*10^12,680*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*680*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*680*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

j = "RCP85outputlowSigma.txt"
dens <- density((UofC - x[2:nrow(x),grep(j,colnames(x))])/scalar)
plot(dens,bty="n",main="",xaxs="i",yaxs="i",xlab="",ylab="",xaxt="n",yaxt="n",type="n",xlim=c(300*10^12,680*10^12)) # xlog=T,xlim=c(10,6*10^14)
polygon(dens, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 20), border=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255, alpha = 200))
abline(v=(UofC - x[1,grep(j,colnames(x))])/scalar, col=rgb(unique(d$r[d$rcp==j]), unique(d$g[d$rcp==j]), unique(d$b[d$rcp==j]), max = 255), )
axis(side=1,at=c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95])), label=round(c(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]),(UofC - x[1,grep(j,colnames(x))])/scalar,max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]))/10^12,0))

# 5% to 95% range is roughly ±5% of mean damages
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar)
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]), y=max(dens$y)/1.3, labels=paste(round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.05]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))
text(x=max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]), y=max(dens$y)/1.3, labels=paste("+", round(100*(max(dens$x[cumsum(dens$y)/sum(dens$y)<0.95]) - (UofC - x[1,grep(j,colnames(x))])/scalar) / ((UofC - x[1,grep(j,colnames(x))])/scalar), 0), "%",sep=""))

# Risk premium
text(x=0.99*680*10^12, y=0.6*max(dens$y), labels=paste("$", round(((x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/10^12,0), "Tr", sep=""), cex=1.2, pos=2)
text(x=0.99*680*10^12, y=0.3*max(dens$y), labels=paste(round((100*(x1[1,grep(j,colnames(x1))] - mean(x1[2:nrow(x1),grep(j,colnames(x1))]))/scalar1)/gdp,0), "%", sep=""), cex=1.2, pos=2)

mtext("Damages (Trillions USD)", side=1, outer=T, line=2.3, cex=.9, at=.5)
mtext(bquote("C=1.2x" * 10^9), side=2,outer=T,las=2,cex=.8, at=.85, line=-1)
mtext(bquote("C=1.0x" * 10^9), side=2,outer=T,las=2,cex=.8, at=.5, line=-1)
mtext(bquote("C=0.8x" * 10^9), side=2,outer=T,las=2,cex=.8, at=.2, line=-1)

mtext(bquote(lambda * "=1.23"), side=2,outer=T,las=2,cex=.8, at=.8, line=.45)
mtext(bquote(lambda * "=1.23"), side=2,outer=T,las=2,cex=.8, at=.45, line=.45)
mtext(bquote(lambda * "=1.23"), side=2,outer=T,las=2,cex=.8, at=.15, line=.45)

mtext(bquote(sigma[Q] * "=1.0625x" * 10^8), side=2,outer=T,las=2,cex=.8, at=.75, line=-2.7)
mtext(bquote(sigma[Q] * "=0.9375x" * 10^8), side=2,outer=T,las=2,cex=.8, at=.4, line=-2.7)
mtext(bquote(sigma[Q] * "=0.6725x" * 10^8), side=2,outer=T,las=2,cex=.8, at=.1, line=-2.7)


dev.off()

